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1.  INTRODUCTION 


,  '  This  paper  studies  two  classes  of  bivariate  regression 
analysis  problems  that  can  be  formulated  and  solved  quite  efficiently 
as  generalized  network  flow  problems. 

Regression  analysis  techniques  are  used  in  many  situations 
to  determine  the  "best"  relationship  between  an  independent  variable 
X  and  a  dependent  variable  Y. 

If  {(xk,yk)}  is  the  set  of  n  paired  observations  of  an 
independent  and  a  dependent  variable,  then  the  bivariate  regression 
problem  is  to  determine  a  linear  regression  equation  of  the  form: 


yk=A+Bxk  (1.1) 

such  that  the  vector  of  predicted  values  of  the  dependent  variable 
is  closa  to  the  vector  of  observed  values.  The  classical  least 
squares  method  of  regression  analysis  uses  the  Euclidian  norm  to 
measure  the  closeness  of  the  vectors  of  predicted  and  observed  values 
[9  ].  Using  this  norm,  the  unweighted  bivariate  least  squares  re¬ 
gression  problem  can  be  stated  as: 

Minimize  (yfc  -  Yk>2^1/2 

subject  to: 

9k  ■  A  +  Bxk  for  k  =  1 ,  2 ,  . . . ,  n. 

This  classical  problem  has  well  known  closed- form  solutions  for  the 
slope,  B,  and  intercept.  A,  of  the  regression  equation  (1.1). 


A 
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The  Euclidian  norm  is  one  member  of  the  L  family  of 

P 

norms.  For  a  specified  value  of  p  =  1,  the  unweighted  bivariate 
1,^  regression  equation  is  given  as  the  solution  to  the  following 
problem : 

Minimize  (yfc  -  yj1^  1/P  (1.2) 

subject  to : 

y  =  A  +  Bx^  for  k  =  1,  2,  . . . ,  n 

The  classical  least  squares,  or  L^,  regression  problem  is 

the  only  member  of  the  L  family  of  regression  problems  that  possesses 

P 

a  closed- form  solution.  In  general,  a  nonlinear  optimization  problem, 

(1.2),  must  be  solved  in  order  to  determine  the  optimal  L  regression 

P 

equation.  Since  this  is  often  an  impractical  means  of  solving  a 
problem,  most  of  the  members  of  the  Lp  family  of  regression  problems, 
other  than  the  classical  least  squares  problem,  have  been  ignored 
in  the  literature.  Two  notable  exceptions  are  the  L^,  or  least  ab¬ 
solute  value,  regression  problem  and  the  L^,  or  Tchebycheff,  regression 
problem. 

Historically,  most  of  the  proposed  solution  procedures  for 
the  (multiple)  regression  problem  have  been  based  on  the  efficient 
simplex  algorithm  of  linear  programming.  However,  the  specific  im¬ 
plementations  of  the  basic  algorithm  have  been  quite  varied.  Chames, 
Cooper,  and  Ferguson  [ 8  )  were  the  first  to  show  that  a  primal  formu¬ 
lation  of  the  problem  could  be  solved  with  the  primal  simplex  algo¬ 
rithm.  Later,  Wagner  [16]  suggested  that  it  would  be  more  efficient 


3 


Co  solve  a  dual  formulation  of  the  problem  since  the  size  of  the 
working  basis  would  be  greatly  reduced.  This  conjecture  was 
widely  accepted  until  Barrodale  and  Roberts  [ 6  ]  showed  that  a 
specialized  primal  simplex  algorithm  could  be  used  to  capitalize 
on  the  underlying  structure  of  the  primal  problem.  McCormick 
and  Sposito  [  11]  proposed  an  improved  Barrodale  and  Roberts  algo¬ 
rithm  which  uses  the  least  squares  regression  problem  to  speed 
convergence.  Abdelmalek  [1,  2]  suggested  that  a  dual  simplex 
algorithm  could  be  used  to  solve  a  dual  formulation  of  the 
problem.  The  equivalence  of  the  Barrodale  and  Roberts  approach 
and  the  Abdelmalek  approach  was  demonstrated  by  Armstrong  and 
Godfrey  [3],  A  specialized  version  of  the  Barrodale  and  Roberts 
procedure  was  developed  by  Armstrong  and  Kung  [5]  to  solve  the 
unweighted  bivariate  regression  problem.  Their  approach  is 

currently  regarded  as  the  most  efficient  one  for  this  class  of 

* 

problems. 

Like  the  problem,  most  of  the  suggested  solution  pro¬ 
cedures  for  the  Lm  (multiple)  regression  problem  have  been  based 
on  the  efficient  simplex  algorithm  [4,  7,  14,  15]. 

This  paper  further  investigates  the  bivariate  and  L^ 
regression  problems.  Specifically,  it  is  shown  that  both  problems 
are  equivalent  to  generalized  network  problems  with  very  special 
network  topologies.  As  a  result,  solution  algorithms  that  exploit 
the  network  topologies  are  developed.  Computer  implementations  of 
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both  algorithms  are  shown  to  be  very  efficient  and  highly  robust. 
This  paper  also  explores  the  underlying  relationship  between  the 
bivariate  and  regression  problems.  It  is  shown  that  the 
problem  is  a  relaxation  of  the  problem,  and  as  such  it  should  be 
easier  to  solve.  This  conjecture  is  supported  by  the  computer  test¬ 
ing  of  the  two  network-based  algorithms. 

2.  WEIGHTED  BIVARIATE  Lj  REGRESSION  ANALYSIS 

Given  a  set  of  paired  observations,  {(x.  and  a 

set  of  n  positive  weights,  {w^},  the  weighted  bivariate  re¬ 
gression  problem  is  to  determine  a  linear  regression  equation 
that  minimizes  the  weighted  sum  of  the  absolute  deviations  (also 
residuals  or  errors).  Formally  stated,  the  problem  is: 

n 

Minimize  I  w.  |  ^  -  y  |  (2.1) 

k-1  K  K 

subject  to: 

y^  *  A  +  Bx^  f or  k  *  1 ,  2 ,  ...,n 

It  is  well  known  that  this  problem  can  be  formulated  as 
a  linear  programming  problem  with  n  +  2  structural  variables  and 
2n  constraints.  This  is  accomplished  by  introducing  a  deviation 
variable,  D^,  and  two  constraints,  -  y^  m  and  ~  yfc)  - 
for  each  of  the  n  pairs  of  observations.  Since  each  weight,  w^,  is 
assumed  to  be  positive,  corresponds  to  the  absolute  deviation 
associated  with  the  kC^  data  point.  That  is, 

\  “  l?k  "  ykl* 
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The  complete  linear  programming  formulation  of  the 
weighted  bivariate  regression  problem  is  given  by: 

n 

Minimize  £  w  D.  (2.2) 

k-1  *  * 


subject  to: 

-A  -  Xj^B  +  1^  m  -  yfc  for  It  -  1,  2 . n  (2.3) 

A  +  XjB  +  S  yfc  for  k  -  1,  2 . a  (2.4) 


All  structural  variables  in  this  primal  linear  programming  problem 
are  treated  as  unrestricted  variables,  even  though  constraints  (2.3) 
and  (2.4)  insure  that  each  deviation  variable,  D^,  is  non-negative. 

The  transformation  of  problem  (2. 2)- (2. 4)  to  its  dual  linear  program¬ 
ming  problem  is  simplified  if  each  deviation  variable  is  considered 
as  an  unrestricted  variable. 

In  order  to  carry  out  the  transformation  of  problem  (2.2)- 
(2.4)  to  its  dual  problem,  it  is  necessary  to  introduce  a  dual  vari¬ 
able  for  each  primal  constraint.  Let  be  the  dual  variable  asso¬ 
ciated  with  the  k1"*1  primal  constraint,  (2.3),  and  let  &k  be  the  dual 
variable  associated  with  the  (n  +  k)^  primal  constraint,  (2.4).  Then 
the  dual  of  problem  (2.2)-(2,4)  is: 

n  n 

Maximize  -  £  yA  +  £  yk&k  (2.5) 

k-1  k-1 

subject  to: 

k-1  “  k-1  “ 

n  n 

"  2  +  2  x  a.  -  o 

k-1  k-1 


n 

I 


(2.6) 


(2.7) 
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\  +  “  Wk  f°r  k  "  1*  2*  •••»  n  (2.8) 

ak  “  0>  8k  =  0  (2.9) 

The  variables  can  be  eliminated  from  the  dual  linear  programming 
problem  by  using  equations  (2.8).  This  is  why  each  of  the  deviation 
variables  in  the  primal  problem  were  treated  as  unrestricted  vari¬ 
ables. 

The  following  problem  is  obtained  from  problem  (2. 5)- (2. 9) 
by  eliminating  the  g^  variables,  scaling  equations  (2.5)  and  (2.7) 
by  -1/2,  and  scaling  equation  (2.6)  by  1/2. 


n 

n 

£  w 

t__1  *■ 

Minimize 

1 

k-1 

K  X 

Vk”  2 

(2.10) 

subject  to: 

n 

n 

-2. 

k-1  * 

(2.11) 

-1 

CM 

1 

R 

k-1 

n 

n 

1 

“  W,  X. 

.  .  ,  k  k 

Vk"  2 

(2.12) 

k-1 

v« 

o 

\-Wk 

(2.13) 

This  formulation  of  the  dual  problem  is  a  capacitated 
generalized  network  problem  with  two  nodes  and  n  arcs.  Figure  1 
illustrates  the  special  network  topology  of  this  problem. 


FIGURE  1 


1^  GENERALIZED  NETWORK 


For  reasons  that  will  become  apparent,  the  node  asso¬ 


ciated  with  constraint  (2.11)  will  be  labeled  A  and  the  node 


associated  with  constraint  (2.12)  will  be  labeled  B.  Node  A 

n 

has  a  net  supply  of  S  w.  /  2  units  and  node  B  has  a  net  demand  of 
n  k=l  k/ 

y  w,x  /  2  units.  It  is  interesting  to  note  that  the  net  gain 
k*l  k  k/ 

through  the  network  is  equal  to  the  weighted  mean  of  observations 


of  the  independent  variable 


k=l 


w  xi 

k  k 


N 


k=l 


All  arcs  for  this  network  problem  are  directed  from  node 
A  to  node  B.  The  kth  arc  has  a  lower  bound  of  zero,  upper  bound 
of  w^,  multiplier  of  x^,  and  objective  function  coefficient  of  y^. 

At  this  point  it  will  be  assumed  that  x.  ^  x.  for  some  i 

i  J 

and  j.  Otherwise,  the  regression  problem  is  meaningless.  This  as¬ 
sumption  is  equivalent  to  assuming  that  the  constraint  matrix  for 
the  generalized  network  problem,  (2. 10)- (2. 13) ,  has  full  row  rank. 

The  generalized  network  problem  always  has  a  finite 
optimal  solution.  This  is  true  since  =  w^y'  2  is  always  a 
feasible  solution,  and  constraints  (2.13)  bound  the  feasible  region. 
Since  problem  (2.10)-(2.13)  always  has  a  finite  optimal  solution, 
duality  theory  ensures  that  the  original  problem,  (2.1),  also 
possesses  a  finite  optimal  solution. 
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Since  problem  (2 .10)- (2. 13)  is  a  capacitated  generalized 
network  problem,  it  can  be  solved  using  a  specialized  network-based 
revised  simplex  algorithm. 

Every  basis  matrix  for  the  generalized  network  problem 
has  the  form: 


B  * 


where  k^  and  k^  are  the  indices  of  the  two  basic  arcs.  The  notation 
B  for  the  basis  matrix  should  not  be  confused  for  the  label  of  the 
node  associated  with  constraint  (2.12).  Since  B  must  have  a  rank 
of  two,  xj^  +  xj^2  *  The  basis  inverse  is 


-1 


'  V 


v 


-k. 


n 


-i 


Let  U  be  an  index  set  of  arcs  that  are  non-basic  at  their 
upper  bounds.  Let 


n 

"I  w 

.  k-1 


I  w, 


keu 


k 


be  the  net  demand  (negative  supply)  at  node  A  after  the  non-basic 
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arcs  have  been  taken  into  consideration.  Similarly,  let 

n 

.  JiV* 

■  '  Jo 


B 


be  the  net  demand  at  node  B. 

The  values  of  the  flows  on  the  two  basic  arcs  are  given 

by : 

f  -  B'-H)  (2.14) 

where 

dk  (\ 

J  V » 

A  network  representation  of  the  basis  is  not  needed  since  (2.14) 
has  the  closed- form  solution: 


bB  +  V  V 


(2.15) 


1  \  ‘  \ 


\  *  \  \ 
C2  ’  \ 


A  basis  is  feasible  if  (2.15)  satisfies  constraints  (2.13) 
for  arcs  k^  and  k^.  The  selection  of  an  initial  feasible  basis  is 
greatly  simplified  by  adding  an  artificial  arc  to  the  generalized 
network  problem.  It  is  convenient  to  number  the  artificial  arc  0. 
Like  the  n  real  arcs,  the  artificial  arc  is  directed  from  node  A  to 


node  B. 
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Let  U  be  the  index  set  of  arcs  such  that  b  <  0.  U  is 

A 

simply  the  set  of  real  arcs  that  are  initially  selected  to  be  non- 
basic  at  their  upper  bounds.  U  can  either  be  initialized  to  be 
empty  or  a  heuristic  rule  can  be  used  to  select  the  arcs  assigned 
to  U. 

Given  an  initial  choice  of  U,  the  multiplier  for  the 
artificial  arc  is  defined  as 


x 

o 


The  artificial  arc  is  given  an  infinite  objective  function  co¬ 
efficient  in  order  to  insure  that  its  flow  in  the  optimal  solution 
is  zero. 

Using  this  artificial  arc,  the  choice  of  an  initial 
feasible  basis  is  very  simple.  Let  k^  *  0  denote  that  the  artifi¬ 
cial  arc  is  the  first  basic  arc.  For  the  second  basic  arc,  any 
arc  k2  t  U  can  be  selected  provided  that  xj^  ^  xj^-  The  flows  on 
these  two  basic  arcs  are  given  by  (2.15).  Due  to  the  specific 

choice  of  x  that  was  made,  the  flows  on  the  two  basic  arcs  are 
o 

simply  ■  -bA  and  ■  0. 

The  specialized  revised  simplex  algorithm  for  solving 
generalized  network  problems  uses  node  potentials  in  order  to  sim¬ 
plify  the  selection  of  entering  arcs.  The  node  potentials  are  de¬ 
fined  as 


ir  -  cbB 


-1 


(2.16) 
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where  u  =  (n^,  t^)  and  =  (y^,  y^)-  A  network  representation 
of  the  basis  is  not  needed  since  (2.16)  has  the  closed- form  solu¬ 
tion  : 


IT 


A 


(2.17) 


7T 


B 


By  applying  duality  theory  to  the  generalized  network 
problem,  (2.10)-(2.13) ,  a  given  basic  feasible  solution  is  an 
optimal  solution  if 


"a  *  V«  ^  \  1£  “k  -  0 
•’a  +  Vb  ;  \  i£  \  •  “k 


If  for  some  arc  k  ,  either 
e 


or 


-IT 


A 


+ 


— TT . 


and 


and 


0 


(2.18) 


(2.19) 


(2.20) 


then  the  current  solution  is  not  optimal  and  arc  kg  is  a  candidate 
to  enter  the  basis  (i.e.,  pivot  eligible). 
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After  selecting  an  arc  to  enter  the  basis  the  specialized 
revised  simplex  algorithm  for  the  generalized  network  problem  must 
compute  the  representation  of  the  entering  arc  in  terms  of  the 
current  basis.  The  representation  of  arc  kg  is  given  by 


The  simplex  algorithm  uses  this  representation  to  perform  the 

standard  minimum  ratio  test.  Specifically,  if  arc  k  is  initially 

e 

non-baaic  with  zero  flow,  then  the  minimum  ratio  is  given  by 


If  arc  ke  is  initially  non-basic  at  its  upper  bound,  then  the  mini¬ 
mum  ratio  is  given  by 


A  ■  min  < 

i  f  \ 

■  w  ,  min<  _Y 

\  <  0 

,  mini 

1 

l  «  i  (  i 

j 

1  i  ( 

If  A  *  w^ ,  then  arc  kg  remains  non-basic  but  it  is  placed  at  its 
opposite  bound.  Otherwise,  arc  kg  becomes  a  basic  arc  and  the  arc 
that  yielded  the  minimum  ratio  is  made  non-basic  at  its  appropriate 
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bound . 

The  status  of  the  variables  in  the  optimal  solution  to 
the  generalized  network  problem,  (2.10)-(2.13) ,  provides  some  use¬ 
ful  insight  into  the  relationship  between  the  original  problem 
data,  {(x^.  y^)},  and  the  optimal  regression  equation.  From  duality 
theory,  it  is  known  that  the  basic  variables  in  the  optimal  solution 
to  the  dual  problem  correspond  to  the  binding  constraints  in  the 
primal  problem.  This  means  that  the  optimal  regression  line  passes 
directly  through  the  two  data  points  corresponding  to  the  optimal 
basis  of  the  generalized  network  problem.  This  result  can  be  used 
to  derive  the  formulas  for  the  slope  and  intercept  of  the  optimal 
regression  equation.  The  slope  of  the  line  through  points 

(Xkx>  and  <*k2’  yk2>  i8: 


B 


(2.21) 


and  the  intercept  of  the  line  is: 


A '  \ '  B\ 


\  yk  "  *k  yk 
*1  *1  *2  *1 

K1  k2 


-IT . 


(2.22) 


It  is  this  relationship  between  the  node  potentials  and  the  optimal 
regression  equation  that  motivated  the  choice  of  A  and  B  as  the  node 
names. 

Using  (2.21)  and  (2.22),  the  optimality  conditions  (2.18), 


become : 
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A  +  Bltk •  yk 

if  Ok  -  Q 

If  \  •  ui 

In  ocher  words,  all  observations  Chat  lie  above  (below)  Che  optimal 
regression  line  have  dual  variables  chat  are  non-baa ic  at  zero 
(upper  bound). 

This  relationship  between  the  sign  of  the  error  term  and 
the  non-basic  status  of  the  arcs  for  problem  (2.10)-(2.13)  can 
be  used  to  develop  heuristics  for  the  initial  selection  of  the  set 
U. 

This  network-based  simplex  algorithm  has  been  coded  in 
FORTRAN  and  tested  on  a  CDC  6600.  Computational  testing  of  this 
code  and  a  state-of-the-art  non-network  code  are  presented  in 
Section  5.  This  testing  indicates  the  superiority  of  the  net¬ 
work  approach. 


3.  BIVARIATE  REGRESSION  ANALYSIS 

Given  a  set  of  n  paired  observations,  {(x^,  y^)),  the 
bivariate  regression  problem  is  to  determine  a  linear  regression 
equation  that  minimizes  the  largest  absolute  deviation.  Formally 
stated,  the  problem  is 


Minimize  max  |  £  -  y 

k  k  k 


subject  to: 


(3.1) 


y^  “  A  +  Bx^  for  k  »  1,  2, 


This  problem  can  be  formulated  as  a  linear  programming 
problem  with  three  structural  variables  and  2n  constraints.  The 


formulation  is 


Minimize 

D 

(3.2) 

subject  to: 

-A  -  xkB 

*  D  -  -*k 

for  k  *  1,  2 . n 

(3.3) 

A  +  XkB 

+  °  -  *k 

for  k  *  1,  2,  . . . ,  n 

(3.4) 

All  three  structural  variables  are  treated  as  unrestricted  variables 
even  though  constraints  (3.3)  and  (3.4)  ensure  that  the  deviation 
variable,  D,  is  non-negative.  Like  the  case,  the  treatment  of 
D  as  an  unrestricted  variable  simplifies  the  transformation  of 
problem  (3. 2)- (3. 4)  to  its  dual  linear  programming  problem.  In 
order  to  carry  out  this  transformation  it  is  necessary  to  introduce 
a  dual  variable  for  each  primal  constraint.  Let  be  the  dual 
variable  associated  with  the  k^  primal  constraint,  (3.3),  and  let 
8^  be  the  dual  variable  associated  with  the  (n  +  k)C^  primal  con¬ 
straint,  (3.4).  The  dual  of  problem  (3.2)-(3.4)  is 

n  n 

Maximize  y  ol  +  I  y  0  v  (3.5) 

k*l  *  *  k*l  *  K 

subject  to: 


n 

I 


n 

+  I 


S. 


(3.65 
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This  new  formulation  is  obtained  by  replacing  equation  (3.6)  with 
one-half  times  equation  (3.6)  less  one-half  times  equation  (3.8), 
replacing  equation  (3.8)  with  one-half  times  equation  (3.6)  plus 
one-half  times  equation  (3.8),  scaling  the  objective  function  (3.5) 
by  -1/2,  and  scaling  each  8^  variable  by  1/x^*  The  scaling  of  each 
3k  implicitly  assumes  that  each  x^  is  positive.  It  will  be  shown 
later  that  this  is  not  a  critical  assumption. 


The  new  formulation  of  the  dual  problem,  (3.10)-(3.14) , 
is  an  uncapacitated  generalized  network  problem  with  three  nodes 
and  2n  arcs.  Figure  2  illustrates  the  special  network  topology 
of  this  problem. 

The  nodes  corresponding  to  constraints  (3.11),  (3.12), 
and  (3.13)  will  be  labeled  A,  B,  and  D,  respectively.  Node  A 
has  a  net  supply  of  1/2  units  and  node  D  has  a  net  demand  of  1/2 
units. 

The  arcs  corresponding  to  the  cx^  variables  are  directed 

from  node  A  to  node  B.  The  multipliers  on  these  arcs  are  x^  and 

the  objective  function  coefficients  are  y^.  The  arcs  corresponding 

to  the  B,  variables  are  directed  from  node  B  to  node  D.  The  multi- 
k 

pliers  on  these  arcs  are  1/x^  and  the  objective  function  coeffi¬ 
cients  are  -y^/x^-  All  2n  arcs  have  lower  bounds  of  zero. 

At  this  point  it  will  be  assumed  that  x,  ^  x.  for  some  i 

k  J 

and  j.  Otherwise,  the  regression  problem  is  meaningless.  This 
assumption  is  equivalent  to  assuming  that  the  constraint  matrix 
for  the  generalized  network  problem,  (3. 10)- (3. 14) ,  has  full  row 
rank. 

The  generalized  network  problem  always  has  a  finite 
optimal  solution.  This  is  true  since  the  network  is  acyclic  and 
*  1/2,  8^  •  x^/2,  and  “  3^  ■  0  for  k  >  1  is  always  a  feasible 
solution. 


) 
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Since  problem  (3. 10)-(3. 14)  is  a  generalized  network 
problem,  it  can  be  solved  using  a  specialized  revised  simplex  algo¬ 
rithm.  Since  the  problem  is  uncapacitated,  at  most  three  arcs  will 
have  non-zero  flows  when  a  simplex  approach  is  used  to  solve  the 
problem.  That  is ,  only  the  three  arcs  associated  with  a  given  basis 
can  have  non-zero  flows.  All  non-basic  arcs  must  have  zero  flow. 

Due  to  the  special  topology  of  the  generalized  network 
problem,  only  two  basis  matrix  structures  are  possible.  For  obvious 
reasons  these  two  structures  will  be  referred  to  as  type  a  and  type 
8  bases. 

The  type  a  basis  consists  of  two  arcs  (k^  and  k^)  and 
one  arc  (k^).  The  type  a  basis  matrix. 


B  - 


0 


-1 


-n 


* 


has  full  row  rank  only  if  xj^  +  x^* 

The  type  8  basis  consists  of  two  6^  arcs  (k^  and  k^)  and 
one  arc  (kj).  cyPe  6  basis  matrix, 


B  - 
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has  full  rank  only  If  xj^  i  x^ • 

Since  B  must  have  full  rank,  It  will  be  assumed  that  the 
basic  arcs  are  numbered  such  that  . 

The  values  of  the  flows  on  the  three  basic  arcs  are 


given  by 


f  -  B-1b 


(3.15) 


where 


-1/2 

b  -  (  0 

1/2 


and  f  “ 


“k! 

^k-2 

ak3 


if  B  is  type  a 


or 


®ki 

ak2 
6k  3  < 


if  B  is  type  B. 


For  a  type  a  basis,  the  solution  to  (3.15)  is  given  by: 


2<\  -  V 


%  '  *k 
1  k2 

2(xk  “  *k  } 
1  3 
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Since  muse  be  non-negative  in  order  for  B  to  be  a  feasible 
basis,  and  since  x^  >  xy^,  it  is  necessary  for  *  xk3* 
Similarly,  since  must  be  non-negative,  »  xj^*  So  B  is 
a  feasible  type  a  basis  if  and  only  if  xj^  a  x^  -  x^  311(1 

xki  >  xk3‘ 

For  a  type  g  basis,  the  solution  to  (3.15)  is  given  by 


\(xk2 '  V 

vv 


\  (xk  “  \  } 

vv 


Since  and  must  be  non-negative,  B  is  a  feasible  type  g 
basis  if  and  only  if  Xj^  «  x^  *  xk3  311(1  xki  >  ^3* 

The  node  potentials  associated  with  a  given  basis  are 


given  by: 

-i 

x  -  cfiB 

where  ir  -  (*A,  ttb,  tfD), 


if  3  is  type  a 


(3.16; 


and 
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or 


if  B  is  type  g. 


For  a  type  a  basis,  the  solution  to  (3.16)  is: 


*A  "  VB  “  \  "  VB  "  \ 

*  VB  *  V 

For  a  type  8  basis,  the  solution  to  (3.16)  is: 


\  •  Vb  "  \ 

nD  "  VB  "  \  ”  V*  "  V 


A  network  representation  of  the  basis  is  not  needed  since  all 
primal  and  dual  variables  have  closed- form  solutions. 

The  value  of  the  objective  function  associated  with  a 


2(xk  '  \  ) 
1  fc3 


The  determination  of  an  initial  feasible  basis  is 
trivial  for  problem  (3.10)-(3.14) .  It  is  only  necessary  to  select 
three  observations,  k][,  k2>  and  k3>  such  that  xk,  >  xk2  l  xk3 
and  xk^  >  xk3<  Let 

9  -  v\  ■  V  ■  V\  •  V +  v\ '  V- 

If  0  ^  0,  let  akl>  gk2  and  ak3  form  the  initial  type  a  basis.  Other¬ 
wise,  let  6^,  ak2  and  Sk3  form  the  initial  type  0  basis. 

Given  a  feasible  basis,  it  is  necessary  to  examine  the 
non-basic  arcs  in  order  to  either  select  an  entering  arc  or  to 
verify  optimality.  From  the  application  of  duality  theory  to  the 
generalized  network  problem,  (3. 10)- (3. 14) ,  a  basis  is  optimal  if 
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ffD  "  Yb  “  yk  ■  *A  for  k  =  1*  2’  •••’  n* 


If  x^TTg  -  then  is  a  candidate  to  enter  the  basis.  If 

x^iTg  -  y^  <  17^,  then  is  a  candidate  to  enter  the  basis.  The 
most  pivot  eligible  arc  can  easily  be  identified  by  solving  the 
following  two  subproblems. 


a 

max  z  ■  x  tt 
k  k  B 


mm  z 
k 


Vb  ' 


(3.17) 

(3.18) 


Ct  S 

If  z  ■  tta  and  zp  ■  ir  ,  then  the  current  solution  is  optimal.  If 

OC  A 

z  -  irA  £  ttd  -  zp,  then  the  a{C£  corresponding  to  the  observation 
that  solved  (3.17)  is  the  most  pivot  eligible  arc.  Otherwise,  the 
Ske  corresponding  to  the  observation  that  solved  (3.18)  is  the  most 
pivot  eligible  arc. 

After  selecting  an  entering  arc,  the  standard  simplex 
minimum  ratio  test  can  be  used  to  determine  the  leaving  arc.  How¬ 
ever,  due  to  the  special  structure  of  the  generalized  network 
problem,  the  minimum  ratio  test  can  be  performed  entirely  with 
logical,  Instead  of  arithmetic,  operations. 

In  order  to  illustrate  how  the  logical  operations  were 
derived,  assume  that  the  current  basis  is  type  a.  If  ot^  is 
selected  to  enter  the  basis,  then  either  or  at^  must  leave  the 
basis  since  a  basis  consisting  of  three  a,  variables  is  not  feas¬ 


ible. 
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Since  xj^  =  xj^  =  xk3»  ak]_  must  leave  the  basis  if 
x^  >  xj^-  This  is  the  only  selection  that  leaves  xj^  in  the 
middle.  Likewise,  if  x^  <  Xj^ ,  then  ot^  must  leave  the  basis. 

If  x^e  =  xk2*  then  either  or  can  leave  the  basis. 

Now  consider  the  case  where  the  initial  basis  is  type  a 
and  3^e  is  selected  as  the  entering  arc.  There  are  three  possi¬ 
bilities  to  consider. 

First,  if  x^e  >  x^,  then  01^3  must  leave  the  basis.  In 
order  to  maintain  xj^  *  ■  xk3’  the  arcs  are  renumbered  as 

follows 


lr  a  Ir 

3  2 

lr  ■  lr 
2  1 

k,  *  k  . 

1  e 

Second,  if  xj^  =  x^  *  xj^ »  then  must  leave  the 

basis . 

Finally,  if  xj^  >  xke>  then  must  leave  the  basis. 
The  arcs  are  renumbered  as  follows 

lr  *  lr 
1  2 

k2  “  k3 

k3  “  ke‘ 


In  the  first  and  third  cases,  the  new  basis  is  type  8. 


In  the  second  case,  the  new  basis  is  still  type  a. 


A  similar  analysis  can  be  done  for  the  case  where  the 


initial  basis  is  type  0. 

After  solving  the  generalized  network  problem,  it  is 
necessary  to  translate  the  optimal  solution  to  the  dual  problem, 
(3. 10)- (3. 14) ,  back  into  an  optimal  solution  for  the  original 
L«)  regression  problem,  (3.1).  By  simply  reversing  the  steps  of 
the  transformation  of  the  dual  linear  programming  problem,  the 
following  intercept  and  slope  of  the  regression  equation  are 
obtained 


A 


♦  V 


B  -  t,b 

The  maximum  absolute  error  is 


As  stated  earlier,  it  is  assumed  that  all  are  positive 
This  is  not  a  critical  assumption  since  the  coordinate  system  can 
be  translated  by  adding  a  constant,  c,  to  each  x^,  such  that  x^  +  c 
is  positive  for  all  observations.  After  solving  the  transformed 
problem,  the  correct  intercept  can  be  obtained  by  adding  c  times 
the  slope  to  the  intercept. 

This  network-based  simplex  algorithm  has  been  coded  in 
FORTRAN  and  tested  on  a  CDC  6600.  Computational  testing  is  pre- 
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senced  in  Section  5. 

4.  RELATIONSHIP  BETWEEN  THE  Li  AND  L*,  DUAL  PROBLEMS 

In  this  section  an  underlying  relationship  between  the 
unweighted  L^  and  L^  regression  problems  will  be  considered.  This 
relationship  is  best  illustrated  by  examining  the  dual  linear  pro-- 
gramming  formulations  of  the  problems. 

The  dual  linear  programming  formulation  of  the  weighted 
L^  regression  problem  was  developed  in  Secion  2.  The  unweighted, 
or  equal  weighted,  L^  regression  problem  can  be  obtained  by  setting 
w^  *  1/n  for  each  observation.  The  unweighted  L^  regression  problem 
is  therefore  equivalent  to 


n 

n 

Maximize 

-I 

k-1 

Vk 

+  ^  yk3k 
k-1  K  K 

(  4.1) 

subject  to: 

n 

r* 

n 

-I 

“k + 

2  6k  -  o 

(  4.2) 

k-l 

k-1 

n 

n 

-I 

k-1 

Vk 

+  1  m 0 
k-1  * 

(  4.3) 

0^+3^“  1/n  for  k  *  1,  2,  ...,  n  (  4.4) 
<^-0,  ■  0  for  k  *  1,  2,  . . . ,  n  (4.5) 

In  Section  3  the  dual  linear  programming  formulation  of 
the  unweighted  L^  regression  problem  was  shown  to  be: 
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Maximize 


subject  to: 


“X  y.  a  +  X  y,8, 

k-1  k-1  K 

n  n 

-X  CL  +  X  8,-0 

k-1  *  k-1  * 


n  n 

-X 


Vk  +  s  *A  “ 0 

k-1  K  K  k-1  K  * 


n  n 

X  \  +  X  3.-1 

k-1  K  k-1  * 


“k  >  0,  ^  >  0  for  k  -  1.  2, 


(4.6) 

(4.7) 

(4.8) 

(4.9) 

(4.10) 


Clearly,  the  dual  problem,  (4.6)-(4.10) ,  is  a  relaxa¬ 
tion  of  the  dual  L^  problem,  (4.1)- (4. 5).  This  relaxation  is 
obtained  by  simply  summing  the  n  constraints  of  form  (4.4).  This 
yields  a  constraint  of  form  (4.9).  Since  the  dual  problem  is 
a  relaxation  of  the  dual  problem,  this  implies  that  the  dual 
problem  should  be  easier  to  solve  than  the  corresponding  dual 
problem.  The  computational  results  presented  in  Section  5 
strongly  support  this  hypothesis  since  the  network-based  code 
runs  one  and  a  half  to  three  times  faster  than  the  network-based 
computer  code. 

5.  COMPUTATIONAL  TESTING 

The  performance  statistics  for  a  network-based  weighted 
regression  computer  code  and  a  network-based  unweighted  re¬ 
gression  computer  code  are  presented  in  this  section.  In  addition. 


for  purposes  of  comparison,  the  solution  times  for  the  state-of- 
the-art,  non-network,  unweighted  regression  computer  code  of 
Armstrong  and  Kung  [5]  are  also  reported. 

All  test  problems  were  generated  by  adding  an  error 
term  to  a  fixed  linear  relationship  between  the  independent  and 
dependent  variables.  The  observations  of  the  independent  variables 
as  well  as  the  error  terms  were  randomly  sampled  from  pre-specified 
uniform  distributions.  Computational  testing,  that  is  not  pre¬ 
sented  here,  indicates  that  none  of  the  solution  algorithms  is 
adversely  affected  by  either  the  range  of  the  independent  vari¬ 
able,  the  magnitude  of  the  error  terms,  or  the  fixed  relationship 
between  the  independent  and  dependent  variables. 

The  non-network  code  of  [5]  uses  a  specialized  dual  simplex 
algorithm  to  solve  a  dual  linear  programming  formulation  of  the  un¬ 
weighted  regression  problem.  The  efficiency  of  this  approach 
is  credited  to  the  fact  that  it  uses  a  technique  referred  to  as 
multiple  pivoting.  Armstrong  and  Kung  have  demonstrated  that  their 
code  is  from  ten  to  one  hundred  times  faster  than  the  code  of 
Sadovski  [12],  which  is  based  on  the  algorithm  of  Edgeworth  [10]. 
Furthermore,  it  has  been  demonstrated  by  Sposito  [13]  that  the 
Sadovski  code  may  fail  to  converge. 

A  discussion  of  the  development  of  the  weighted  re¬ 
gression  code  will  be  presented  next.  This  is  following  by  a  compu¬ 
tational  comparison  of  this  network-based  code  to  the  state-of-the- 
art  code  of  Armstrong  and  Kung.  Then  a  discussion  is  given  on  the 
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impact  of  the  magnitude  and  range  of  the  weights  on  the  weighted 
regression  code.  After  presenting  the  code,  the  development 
and  computational  testing  of  the  regression  code  will  be  out¬ 
lined.  This  is  followed  by  a  comparison  of  the  two  new  network 
codes . 

Two  of  the  most  crucial  determinants  of  the  relative  ef¬ 
ficiency  of  network-based  simplex  codes  are  (1)  the  choice  of  the 
starting  basis  (i.e.,  the  rule  for  constructing  the  initial  basic 
feasible  solution)  and  (2)  the  choice  of  the  pivot  strategy  (i.e., 
the  rule  for  selecting  entering  arcs).  As  expected,  both  of  these 
aspects  have  major  ramifications  for  the  solution  efficiency  of 
the  weighted  regression  code.  Two  start  rules  were  developed 
and  tested  as  well  as  two  standard  pivot  rules.  A  brief  discussion 
of  each  follows. 

In  Section  2  It  was  shown  that  an  initial  feasible 
(artificial)  basis  can  be  constructed  by  (1)  selecting  a  set  of 
arcs  U  to  be  non-basic  at  their  upper  bounds,  (2)  constructing  a 
basic  artificial  arc  that  transfers  the  remaining  units  of  supply 
from  node  A  to  node  B,  and  (3)  selecting  a  real  arc  to  complete 
the  full  row  rank  basis.  The  most  important  step  of  this  start 
procedure  is  the  Initialization  of  the  set  U.  The  first  start  rule 
tested  simply  initializes  U  as  the  empty  set.  The  second  start 
rule  that  was  tested  was  motivated  by  an  observation  based  on 
duality  theory.  This  observation  is  that  ^  £  y^  (^  ■  y^)  if  the 
optimal  network  solution  has  oi^  ■  0  (a^  -  w^).  The  heuristic  start 


rule  based  on  this  observation  is  referred  to  as  the  advanced 
start.  Quite  simply,  the  heuristic  is  to  use  the  weighted  least 
squares  regression  equation  to  assign  arcs  to  the  set  U.  That  is, 
U  initially  contains  those  arcs  associated  with  positive  least 
squares  residuals.  The  solution  times  (in  c.p.u.  seconds)  and  the 
required  pivots  for  both  start  rules  are  presented  in  Table  1. 

All  testing  was  done  on  The  University  of  Texas'  CDC  6600  using 
the  MNF  FORTRAN  compiler.  In  order  to  reduce  the  effect  of 


TABLE  1 
START  RULE 


EMPTY  U 

ADVANCED  L2 

NUMBER  OF 
OBSERVATIONS 

SOLUTION 

TIME 

PIVOTS 

SOLUTION 

TIME 

PIVOTS 

100 

.174 

102 

.016 

8 

200 

.671 

203 

.049 

11 

300 

1.474 

301 

.148 

29 

400 

2.604 

400 

.166 

23 

500 

4.025 

498 

.291 

32 

*L0NE  with  most  eligible  pivot  rule 
Mean  times  and  pivots  reported 


measurement  error,  all  solution  times  and  pivots  that  are  reported 
are  actually  averages  for  numerous  test  problems. 

It  is  clear  from  Table  1  that  the  advanced  start  is 
far  superior  to  the  simple  empty  U  start.  It  is  important  to 
note  that  the  time  required  to  determine  the  weighted  least  squares 
regression  equation  is  included  in  the  reported  times  in  Table  1 
for  the  advanced  start  code,  but  is  not  included  for  the  empty 
U  start  since  it  is  not  used  by  the  code. 

The  second  aspect  of  Che  network.- based  regression  code 
that  was  studied  was  the  choice  of  pivot  rule.  Numerous  rules  have 
been  suggested  in  the  network  literature  for  the  selection  of  enter¬ 
ing  arcs.  Two  of  the  most  straightforward  rules  are  to  pivot  on 
(1)  the  most  pivot  eligible  arc,  and  (2)  the  first  pivot  eligible 
arc  found.  The  philosophy  underlying  these  rules  is  that  the  most 
eligible  rule  requires  a  lot  of  effort  to  identify  the  entering 
arc,  but  is  expected  to  produce  the  fewest  total  pivots.  On  the 
other  hand,  the  first  eligible  rule  requires  little  effort  to 
identify  an  entering  arc,  but  is  expected  to  require  more  total 
pivots.  Both  pivot  strategies  were  tested  in  order  to  determine 
the  relative  trade-off  between  the  amount  of  work  required  and 
the  number  of  pivots  performed.  Table  2  presents  the  solution 
times  and  pivots  for  these  rules.  Again,  the  reported  statistics 
reflect  the  averages  of  numerous  test  problems. 


The  impact  of  the  pivot  rule  is  clear.  The  first  eli¬ 
gible  rule  requires  from  two  to  four  times  as  many  pivots  as  the 


PIVOT  STRATEGY  RULE 


MOST  PIVOT 

ELIGIBLE 

FIRST  PIVOT  ELIGIBLE 

NUMBER  OF 

SOLUTION 

SOLUTION 

OBSERVATIONS 

TIME 

PIVOTS 

TIME 

PIVOTS 

50 

.016 

11 

28 

100 

.016 

8 

.016 

16 

150 

.032 

9 

31 

200 

.049 

11 

-  BS m 

33 

250 

.114 

24 

60 

300 

.148 

29 

.063 

72 

350 

.179 

28 

.065 

61 

400 

.166 

23 

.088 

88 

450 

.300 

37 

.102 

107 

500 

.291 

32 

.108 

123 

l  -  -  - ; 

*LONE  with  L2  advanced  start 
Mean  times  and  pivots  reported 
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mosc  eligible  rule.  However,  the  first  eligible  rule  is  consis¬ 
tently  faster  than  the  mosc  eligible  rule. 

Based  or.  the  computational  testing  of  the  start  rules 
and  the  pivot  rules,  the  be't  network-based  weighted  regression 
code,  referred  to  as  LONE,  uses  the  advanced  start  and  the  first 
eligible  pivot  rule.  In  order  to  properly  assess  the  performance 
of  this  network  approach,  the  fastest  non-network  L^  code  known 
was  acquired.  This  code,  known  as  SIMLP,  is  based  on  the  dual 
simplex  algorithm.  Unlike  LONE,  SIMLP  can  only  solve  unweighted  L^ 
regression  problems.  Both  codes  require  four  n  length  data  arrays, 
but  a  streamlined  unweighted  version  of  LONE  could  be  developed 
chat  uses  only  three  n  length  data  arrays.  The  solution  speed  of 
such  a  code  would  also  improve  since  the  operations  using  the 
weights  would  be  streamlined. 

Table  3  presents  the  computational  results  for  LONE  and 
SIMLP  on  a  large  set  of  unweighted  test  problems.  The  indication 
is  that  LONE  is  roughly  twice  as  fast  as  the  state-of-the-art  code 
SIMLP.  This  supports  the  recent  findings  that  network  simplex  algo¬ 
rithms  are  computationally  superior  to  their  non-network  counter¬ 
parts. 

Limited  computational  testing  of  LONE  was  carried  out  in 
order  to  measure  the  impact  on  its  performance  capabilities  as  the 
magnitude  and  range  of  the  weights  are  varied.  Table  4  reports 
the  average  solution  times  and  pivots  for  a  number  of  test  prob¬ 
lems  with  300  data  observations.  For  each  of  these  problems,  the 


TABLE  3 

LONE  VS.  SIMLP 


LONE 

SIMLP 

NUMBER  OF 

SOLUTION 

SOLUTION 

MULTIPLE 

OBSERVATIONS 

TIME 

_ 

PIVOTS 

TIME 

PIVOTS 

50 

mm 

28 

.010 

3 

100 

16 

.023 

4 

150 

■l/ 

31 

.038 

6 

200 

MEM 

33 

.049 

4 

250 

■Si' 

60 

.076 

4 

300 

.063 

72 

.095 

5 

350 

.065 

61 

.129 

5 

400 

.088 

88 

.139 

4 

450 

.102 

107 

.180 

5 

500 

.108 

123 

_ 

.213 

6 

*LONE  with  L2  advanced  start  and  first  eligible  pivot  rule 
Mean  times  and  pivots  reported 


TABLE  4 

WEIGHTED  PROBLEMS 


MINIMUM 

WEIGHT 

MAXIMUM 

WEIGHT 

SOLUTION 

TIME 

PIVOTS 

1 

1 

.064 

76 

1 

2 

.064 

70 

1 

5 

.059 

74 

1 

10 

.057 

63 

*300  observations 
Mean  times  and  pivots  reported 
LONE  with  L^  start  and  first  eligible  pivot 
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weights  were  uniformly  distributed  between  a  specified  upper  and 
lower  bound.  This  testing  indicates  that  LONE  is  insensitive  to 
the  variation  in  the  weights. 

A  brief  summary  of  the  development  and  computational 
testing  of  the  network-based  unweighted  L^  regression  code  is  pre¬ 
sented  next.  As  is  the  case  with  all  network  simplex  codes,  the 
most  critical  algorithmic  decisions  made  during  the  development 
of  the  L  code  dealt  with  the  choice  of  the  initial  basis  and  the 

OQ 

entering  arc  selection  rule.  The  most  efficient  implementation, 
referred  to  as  LINF,  constructs  its  initial  basis  from  the  first 
three  observations,  (xkl,  ykl>,  (xk2>  yk2>,  and  (xk3>  yk;}) ,  such 
that  xkl  +  xk;J.  This  code  uses  the  most  eligible  rule  to  select 
the  entering  arcs.  LINF  only  uses  two  n  length  data  arrays. 

The  solution  times  and  the  number  of  pivots  performed  by 
LINF  are  reported  for  a  large  set  of  unweighted  test  problems  in 
Table  5.  For  purposes  of  comparison,  the  table  also  provides  the 
results  for  the  network-based  code  on  the  same  problems.  As 
hypothesized  in  Section  4,  the  L^  problem  is  indeed  easier  to 
solve  than  the  corresponding  L^  problem.  This  is  due  to  the 
fact  that  the  L^  dual  linear  programming  problem  can  be  inter¬ 
preted  as  a  relaxation  of  the  L^  dual  linear  programming  problem. 
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TABLE  5 

LINF  VS.  LONE 


LINF 

LONE 

NUMBER  OF 
OBSERVATIONS 

SOLUTION 

TIME 

PIVOTS 

SOLUTION 

TIME 

PIVOTS 

50 

.004 

3 

.016 

28 

100 

.008 

4 

.016 

16 

150 

.010 

5 

.015 

31 

200 

.015 

4 

.022 

33 

250 

.026 

5 

.053 

60 

300 

.022 

5 

.063 

72 

350 

.035 

6 

.065 

61 

400 

.050 

6 

.088 

88 

450 

.046 

6 

.102 

107 

500 

.043 

6 

.108 

123 
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